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Abstract Realistic networks display heterogeneous transmission delays. We analyze here the limits 
of large stochastic multi-populations networks with stochastic coupling and random interconnection 
delays. We show that depending on the nature of the delays distributions, a quenched or averaged 
propagation of chaos takes place in these networks, and that the network equations converge towards a 
delayed McKean-Vlasov equation with distributed delays. Our approach is mostly fitted to neuroscience 
applications. We instantiate in particular a classical neuronal model, the Wilson and Cowan system, 
and show that the obtained limit equations have Gaussian solutions whose mean and standard deviation 
satisfy a closed set of coupled delay differential equations in which the distribution of delays and the 
noise levels appear as parameters. This allows to uncover precisely the effects of noise, delays and 
coupling on the dynamics of such heterogeneous networks, in particular their role in the emergence of 
synchronized oscillations. We show in several examples that not only the averaged delay, but also the 
dispersion, govern the dynamics of such networks. 

Keywords Heterogeneous neuronal networks • Mean-field limits • Delay differential equations • 
Bifurcations 
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Introduction 

Realistic networks show heterogeneous interconnection delays. The cortex is a paradigmatic example of 
such networks. In the brain, neurons form large-scale assemblies (termed neural populations) subject 
to an intense noise. They communicate through the emission of spikes (action potentials) that are 
transmitted through synapses that can be chemical or, less often, electrical. The transport of spikes 
through the axons takes a specific time, function of the length of the axon and of the time taken for the 
signal to be transmitted through synapses, and therefore are highly heterogeneous. Moreover, chemical 
synapses modify the membrane potential of the post-synaptic neuron through a complex mechanism of 
release and binding of neurotransmitters. This mechanism takes a specific time, and results in random 
variations of the current transmitted through the different synapses and also from spike to spike. 

Models describing the emergent behavior arising from the interaction of neurons in large-scale 
networks largely overlooked both phenomena of heterogeneous delays and stochastic synapses. Most 
models have relied on continuum limits ever since the seminal work of Wilson and Cowan and Amari 
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2.3, 32 , 33J . Such models tend to represent the activity of the network through a macroscopic variable, 
the population-averaged firing rate, that is generally assumed to be deterministic. Many analytical 
and numerical results and properties have been derived from these equations and related to cortical 
phenomena, for instance for the problem of spatio-temporal pattern formation in spatially extended 
models (see e.g. [5l [T2in51[T7] ) . This approach tends to implicitly make the assumption that the effect 
of noise and heterogeneities vanishes in large populations. A different approach has been to study 
regimes where the activity is uncorrelated. A number of computational studies on the integrate-and-fire 
neuron showed that under certain conditions neurons in large assemblies end up firing asynchronously, 
producing null correlations [1,4.5 . In these regimes, the correlations in the firing activity decrease 
towards zero in the limit where the number of neurons tends to infinity. The emergent global activity 
of the population in this limit is deterministic, and evolves according to a mean-field firing rate equation. 
However, these states arise in sparsely connected networks which is not the case of cortical columns. 
In order to go beyond asynchronous states and take into account the stochastic nature of the firing 
and how this activity scales as the network size increases, different approaches have been developed, 
such as the population density method and related approaches [7j. Most of these approaches involve 
expansions in terms of the moments of the resulting random variable, and the moment hierarchy 
needs to be truncated which is a quite hard task that can raise a number of technical issues (see 
e.g. [18]). Homogeneous spatially extended systems of neurons were recently analyzed using the theory 
of probability [30ll29| . These were the first studies considering individual delays between neurons. In 
these studies, the propagation delays between neurons were considered constant, and in these cases, no 
specific delay-induced transition was identified. Heterogeneous networks have been considerably less 
analyzed. One notable exception is the work of Sompolinsky and collaborators |25j recently analyzed 
from a mathematical viewpoint using large deviations techniques [6], where the authors considered 
random heterogeneous connections between neurons. 

The main question we shall address in this article is the effect of heterogeneous delays on the 
macroscopic dynamics of such neuronal areas in the presence of stochastic coupling. Although it now 
firmly established that delays are fundamental for shaping the activity of large-scale neuronal net- 
works 22.8,23,24], very little is understood about the type of delays involved in such networks. The 
question of how these individual delays translate in macroscopic limits, and how the heterogeneity in 
the distribution of delays affect macroscopic behaviors, hence has deep implications, and is still largely 
open. At the microscopic scale, it is relatively clear that transmission delays are well approximated 
by constant delay depending on the distance between neurons and on the typical synaptic time scale. 
These constant delays hence vary depending on the pair of neurons considered and are related to the 
topology of the network in the physical space. 

Considering random heterogeneous delays between neurons, we establish the limit, as the number 
of neurons tend to infinity, of the network equations. We will show that the propagation of chaos 
property occurs in these systems, and that the asymptotic macroscopic behavior in the limit where the 
number of neurons tends to infinity is described by a McKean-Vlasov equation with distributed delays, 
the delay measure being exactly the statistical distribution of individual delays. The dynamics of the 
obtained McKean-Vlasov equations with distributed delays is very hard to analyze for general neuron 
models. In order to uncover the macroscopic dynamics of the neurons and their dependence on the 
distribution of delays, we eventually instantiate a particularly suitable model, the classical firing-rate 
model. In that case, that solutions are Gaussian, and the mean and standard deviation satisfy a closed 
system of delay differential equations with distributed delays. This reduction allows using classical 
tools of the domain of infinite-dimensional dynamical systems to show generic qualitative transitions 
between different qualitative states, as a function of the delays and of the variance of the synaptic 
weights. 

1 Setting of the problem 

In this section, we describe the mathematical formalism used throughout the manuscript. We shall 
analyze the dynamics of a relatively general neuron model, valid for most usual models used in com- 
putational neuroscience. The state of each neuron i in the network is classically described by a d- 
dimensional variable X 1 EE:— R d , typically corresponding to the membrane potential of the neuron 
and possibly additional variables such as those related to ionic concentrations and gated channels (see 

e.g. m). 
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We consider networks composed of a finite number of populations a G {1, • • ■ , P}. The populations 
differ through the intrinsic properties of each neurons and the input they receive. We assume that 
the number of neuron in each population a G {1, . . . , P}, denoted N a increases as the network size 
increases, and tends to infinity in the limit considered, which we will be denoting N — > oo. We define 
the population function that associates to a neuron i the population a it belongs to: p(i) = a. 

Let us consider a single neuron i belonging to population a. The state of the neuron, X 1 , has an 
intrinsic dynamics governed by a nonlinear drift function f a : R x E t— > E (including the intrinsic 
dynamics and the deterministic inputs) and a diffusion matrix g a :Rx£4 E m . Each individual 
process {X l t ) takes values in E, and when there is no interaction, is governed by the equation: 

dX\ = f a (t, XI) dt + g a (t, XI) dWl 

where (W t *)j g ]N are independent m x d Brownian motions. 

This equation corresponds to the intrinsic dynamics of each individual diffusion, i.e. the behavior of 
a single neuron when disconnected of the network. When included in the network, the neurons interact 
with all other neurons and this interaction is modeled through a continuous function b ai {x, y) : ExE i-» 
E, multiplied by a coefficient, called the synaptic weight w ir These synaptic weights are themselves 
stochastic, due to the nature of the synaptic transmission. Incoming interactions are known to be 
bounded, the total incoming connectivity from a given population (say 7) to a single neuron from 
population a has a fixed averaged J ai and standard deviation a ai . We consider that the synaptic 
weights are stochastic, and are assumed to be equal, for a neuron i in population a and a neuron j in 
population 7, to 

w=w =—J + — B n 

where B 11 are independent <5-dimensional "white noise" processes related to the independent Brownian 
motions B %a introduced^! 

The delay r,-j occurring in the interaction between neuron j and i, and related to the information 
transmission through the axons and the synapses, depend on the precise distance between neurons i 
and j. We assume that for any fixed neurons i G {1 ■ • • N}, the delays (t^)^ are pairwise indepen- 
dent (and independent of the Brownian motions involved in the network dynamics) and identically 
distributed in each population, with a law denoted This is a reasonable assumption from the 

modeling viewpoint. Indeed, neural populations can be considered to form certain clusters in the space 
in which the specific locations of neurons can be considered independent, and delays hence depend on 
the typical distance between the populations, their dispersion, and on the characteristic time constants 
of the synapse. From the referential of one given neuron i, the resulting delays are hence independent. 
From the global network viewpoint, these delays are of course correlated (distance between two points 
are symmetrical, topological relationship between neuronal locations induces relationship between dis- 
tances, e.g. triangular inequality). Let us denote by r the largest possible delay in the network. All 
delay distributions hence charge only the interval [— r, 0]. We eventually denote by r\ ai the averaged 
delay between neurons of population a and neurons of population 7, were the average is taken over all 
possible values of rj^ for i in population a (for instance, averaged on all possible locations for neurons 
of population a). 

The infinitesimal current the neuron i receives from its neighbors is hence modeled as: 



7 =1 p(j)=l 



N. 



P(j)=7 



{X, 



i.N 



X 



sdBV 



1 Note that this model has the disadvantage to possibly change the excitatory or inhibitory nature of the 
connection from neuron i to neuron j when the sign of Wy changes. Other more biologically relevant problems 
involve enforcing the synaptic weight not to change sign by modeling it as the solution of particular stochastic 
differential equation that does not change sign, such as the Cox-Ingersoll-Ross process, and by fit our framework 
by formally adding the P synaptic weights of neuron i in the state variable X 1 of neuron i. 
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For the sake of simplicity, we make an abuse of notations and replace the notation J ai b ai {x,y) 
(resp. a ai j3 ai {x,y)) by b a7 (x,y) (resp. j3 aj (x,y)). The resulting network equation reads: 

p N 

dXt> N = f a (t,Xi< N )dt + J2 E ^b ai {X[- N ,Xt N Tij )dt 

7=1 j=l,2>0')=7 7 

P AT 

+ ga (t,xi> N )-dwt+j2 w E /m^^)-^ 7 - a) 

7=1 7 J=l,p(j)=7 

These equations are stochastic differential equations on the infinite-dimensional space of functions 
C([— t, 0], E) (i.e. on the variable X t = (X s ,s £ [t — r,i\), see e.g. [1THI20] ) of continuous functions of 
[— r, 0] with values in E. Let us denote by C T = C([—r, 0], E p ). The initial conditions are not specified 
at this point. We will sometimes make the assumption that the network has chaotic initial condition, 
i.e. independent identically distributed initial conditions. In details, for (C^(t),o; = 1 • • ■ P) £ C T a 
stochastic process with independent components, a chaotic initial condition on the network consists in 
setting the initial condition of all neurons i of population a to an independent copy of Co*- 

For any (a, 7) € {1, • ■ ■ , P} 2 , we make the following technical assumptions on the parameters of 
the model: 

(HI). f a and g a are uniformly locally Lipschitz-continuous with respect to their second variable 
(H2). b aj and /3 ai are L-Lipschitz-continuous, i.e. there exists a positive constant L such that for all 
(x, y) and (2', y') in E x E we have: 

\0 ai (x, y) - O ai (x', y')\<L (\x - x'\ + \y- y'\) 

for O € {b, /?}. 
(H3). There exists a K > such that: 

max(|6 Q7 (a;,z)| 2 , \ j3 ai (x, z)\ 2 ) < K 
(H4). The drift and diffusion functions satisfy the monotone growth condition: 

x T f a (t,x) + l -\g a {t,x)\ 2 <K(l + \x\ 2 ). 

Remark 1 Assumptions [(HI)] and [(H4) | are technical assumptions necessary to deal with the particular 
case of the Fitzhugh-Nagumo neuron model [15]. Most classical neuron models are however defined 
through globally Lipschitz-continuous vector fields. When the drift and diffusion functions are only 
locally Lipschitz continuous, a classical localization argument allows extending the proofs done in the 
globally Lipchitz continuous case, as done in |30j . This is why in the manuscript we will consider the 
case of globally Lipschitz continuous functions f a and g a and refer the interested reader to |30j for 
extending these proofs to locally Lipschitz continuous drift and diffusion functions. 

Let us first state the following proposition ensuring well-posedness of the network system under 
the assumptions of the section: 

Proposition 1 Let Xq £ Ai 2 (C([— t, 0], E N )) an initial condition of the network system. Under the 
assumptions of the section, there exists a unique strong solution to the network equations ([T]). 

The proof of this proposition is a direct application of Mao's result [T5], and essentially uses the 
same arguments as those of the proof theorem [T] The interested reader is invited to follow the steps 
of this demonstration to prove proposition [T] 

We address the problem of the asymptotic behavior, as the number of neurons N goes to infinity, 
of these network equations, in particular the question of the propagation of chaos property in this 
randomly delayed, randomly interacting system. Assuming translation invariance in the neural field, 
the delays distribution rjij is only function of the neural population i belongs to. In that case, for any 
realization of the synaptic delays ry, the propagation of chaos property states that, provided that the 
initial condition on the network is chaotic (i.e. for all i in population a, the initial condition JQ|[_ r ] 
is an independent copy of a given process € C T ) and when N — > 00, the state of each neurons 
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are independent for all times, and have a law given by a nonlinear McKean-Vlasov equation that 
only depends on the neural population (say, a), which is the unique solution of the McKean-Vlasov 
equations with distributed delays: 



dX? = f a (t,X?)dt + J2[ E y [b afi [X? , Yf +s )]d Va p (s)dt 

P M 

+ g a (t, X«) dW? + / ^Y[f3 afS (X^Yf +s )]dvMs) dBf (2) 

7=1 J - T 

where Y is a process independent of X that has the same law, Ey the expectation under the law of 
Y, W" and B" 7 are independent adapted standard Brownian motions of dimensions m x d and S x d 
respectively. In other words, the mean field equation can be written, denoting by mj(dx) the law of 

X7-. 



y)mj +s (dy) dr] aj (s)dt 



dX? = f a (t, X? )dt + J2 J «l / / bay (X? , 
7=1 J-tJE 

+ 9a (t,X?)dW t a + J2j ai f I 

7=1 J-tJE 



Pa-fiX?, y)m] +s (dy) drj ai (s)dB? 



where m" is the probability distribution of X". In these equations, W° for a = 1 • • • P are independent 
standard m-dimensional Brownian motions. This equation is also classically written as the McKean- 
Vlasov Fokker-Planck equation on the probability density, with a slightly abuse of notations denoted 
m](y), of the solution: 



d t m?(x) = -V^/ a (t,i) m?(xf) -^x(^2j j b ai (x,y)m] +s (y) dy dr] a7 (s)m?(x)^j 

+ -A x (g a {t,xfm?{x)) + -A x { ^ ( j J J ai (x,y)m< +S {y) dy dry Q7 ( S )) m?(x)} (3) 



7=1 



with initial condition mf(x)\t£[^ r ,o] = P a (t, x) where p a is the density of Co 7 if it exists. 

The same result holds when the neural field is not translation invariant (i.e. 77^ varies depending 
on the choice of neuron i in population a), but in that case the result holds for the law of the limit 
distribution averaged on the law of the delays. 

In that view, the equation ^ is an implicit equation on the law of X t (denoting when no confusion 
is possible the vector (X t Q , a = 1 • • ■ P)). In order to prove the propagation of chaos property, we will 
use the coupling method (see e.g. [27]). The proof is in two steps: (i) we prove that the equation ^ has 
a unique solution, and (ii) that the law of X l t ' N converges towards the law of |2| and more precisely 
that for any finite set of neurons {it, ■ ■ ■ ,ik} the law of (X t ll,Ar , • • • ,X\ h ' ,t G [— t,T]) converge almost 
surely towards a vector (Xh--- ,X^,t G [— r, T]) where the processes X 1 are independent and have 
the law of given by §. 



2 Well-Posedness of the delayed McKean-Vlasov equations 

In this section, we show existence and uniqueness of solutions for the delayed McKean-Vlasov equa- 
tion (2§ that constitute the putative mean-field limit. Let us denote by M(C) the set of probability 
distributions on C the set continuous functions [— r, T] E p , and A4 2 (C) the space of square- integrable 
processes. 

Let (W a ; a = 1 • • • P) (respectively (B" 1 , ; a, 7 = 1 • • • P)) be a family of P (resp. P 2 ) inde- 
pendent m x d (resp. 6 x d)-dimensional standard Brownian motion adapted on (Q,F,P). Let us also 
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fixe Co € Ai 2 (C T ) a random variable, and (Co)ign an iid sequence with the same law as Co' 1 '*' that 
constitutes the chaotic initial condition of the network equations. We introduce the map <P map acting 
on stochastic processes and defined by: 

'M{C) ^M(C) 

X ^(Y t = {Y t a ,a = l--- P}) t with 

y t a = Co a (0) + Jo (Us, X?) + Ej =1 fir Ez[6 Q7 (X«, Z2 +u )]) di lai {u) ds 
+ Si g a (s, Xf) dW« + Ef =1 So Sl T Z] +u )} d Vai (u) ■ dBf ; t > 



Y t = Q(t) t€[-r,0] 

where the process introduce a process Z t has the same law as X and to be independent of X. There 
is a trivial identification between the solutions of the mean-field equation ([2| and the fixed points of 
the map <£: any fixed point of <P provides a solution for equation ^ and conversely any solution of 
equation ([2| constitute a fixed point of <£. 

Let us start by proving that the possible solutions of the system have bounded second moment. 

Lemma 1 Let € Ai 2 (C([— t, 0], E p ) a square-integrable process. Under the assumptions (H1)-(H4), 
there exists a constant C{T) > depending on the parameters of the system and on the horizon T, 
such that for any solution X of the mean-field equation |2| with initial condition (q: 

sup n\Xt\ 2 } <C(T). 

l-r,T] 

Proof The proof is based on the application of Ito's formula for the squared modulus of X tl standard 
inequalities and Gronwall's lemma. 

Let X be a solution of the mean-field equations, and r„ the first time the process \X t \ exceeds 
the quantity n. Due to the non-standard nature of the equation, let us underline the fact that Ito's 
formula is valid, i.e. that X t is a semimartingale. By definition, it is clear that both X t + S and Z t + S 
are !F t measurable for all s£ [— r, 0], and hence necessarily X t is a semi-martingale, i.e. the sum of a 
continuous adapted process of finite variation: 

jf (/«(*, X a s ) + f 1MW*7> Z] +u )] dr> aj (u)) ds 

and a continuous (P t , P)-local martingale, defined as the sum of different stochastic integral of a 
progressively measurable processes with respect to Brownian motions: 

f g a {s,X*)dW« + Yl f f Vz[(3 a7 (Xs,Z] +u )]d Va7 (u)-dBr 

JO 7 =l^ -'- T 

We can hence apply ltd formula to |X tATi J 2 , and because of the particular form of the squared norm, 
we can treat each component a £ {!,••• , P} separately. Since all the different Brownian motions 
involved are independent, we have for all t > 0: 



\X^rf 



IC Q (o)l 2 



(I 



{(X«) T f a (s,X?)+'l\g a (s,X?)\ 2 



+ PC) T ]T f M Y {b a ^,Y s \ u )]d Vaj (u) + ^\ f ^Y[(3 a y(Xf,Y s l u )]drj a7 (u)\ 2 } 

7=1 J - T 7=1 J - T 



ds 



t/\T n 



\x«) T g a (s,X«)dW t a + J2 



7=1 



*At„ 



m Y \p an {xf,Y2^)\d ncei {u)dB_ 



The stochastic integral term has a null expectation, and the Stie fjes in tegration term involves the 
term \x T b a ^(x 7 z)\, which is upperbounded because of assumption (H3) by + \x\ 2 ), and the 
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term x T f a (t, x) + h \g a (t, x)\ 2 which is upperbounded, using assumption (H4) by K(l + \x\ 2 ). Finally, 



assumption 
hence have 



(H3) | again allows us to upperbound the term ±E y [\(3 ai {Xf , Y?)\ 2 } by f (1 + \Xf\ 2 ). We 
or all t > 0: 



E[|X« T J 2 ]<E[|C o a (0)| 2 ]+2 J {X(1+E[|^| 2 ]) + PVX(1+E[|X S Q | 2 ]) 

+ yE f° a + n\K\ 2 ])dv a ,(u)}d s 

7=1 J ~ T 



<E[|C Q (0)| 2 ]+2 J (K+^KP+-P)(\+n\Xs\ 2 ])ds 



Summing over a and applying Gronwall's lemma directly yields: 



sup E[|X t | 2 ] < E[|Co(0)| 2 ] +E[1 + |Co(0)| 2 ](e K ' T - 1) 

tG[0,TAr„] 



with K' = 2 [K + #P + V K P), and hence, letting n — > oo, we conclude that: 



sup E[|X t | 2 ] <max(sup E[|Co(s)| 2 ] , E[|Co(0)| 2 ] + E[l + |Co(0)| 2 ](e K ' T - 1)) 
te[-r,T] [-r,0] 



We now prove the existence and uniqueness theorem, using this boundedness property of the pos- 
sible solutions. 



Theorem 1 For any Co G A4 2 (C([—t,0], E p ) a square-integrable process, the mean-field equation ([2| 
with initial condition Co has a unique strong solution on [— r, T] for any T > 0. 



Proof We start by showing the existence of solutions, and then prove the uniqueness property. We recall 
that by application of lemma[TJ any possible solution has a bounded second moment. The proof is based 
on a classical contraction argument on iterates of the map There are two main differences with usual 
proofs. First is the fact that our functions are not globally Lipschitz-continuous, issue which is solved 
through a truncation argument and using the result of the a priori estimate of lemma [lj see |30j . 
Second is the fact that we are considering delayed equations, which leads us to deal with uniform in 
[— r, T] convergence properties. The rest of the proof uses classical inequalities such as Cauchy-Schwarz 
(CS), Burkholder-Davis-Gundy (BDG) and Holder inequalities, all these applied on a decomposition 
of the term of interest into elementary terms. 
Existence: 

Let X° E Ai 2 (C) such that X°|[_ t ] =Co a given stochastic process. We consider that the drift and 
diffusion functions are ii"-Lipschitz continuous (up to truncation). 

We introduce the sequence of probability distributions (X k )k> on Ai(C) defined by induction 
as X k+1 = (<P(X k )). We denote by (Z k ) a sequence of processes independent of the collection of 
processes (X k ) and having the same law. We start by controlling one of the components of the sequence 
of processes (X k ), and for compactness of notations denote a X k e E the component a of X k . We 
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decompose into six elementary terms the difference: 

a x k+l _ aj{ k = J* ^ /q(S; a x k } _ a X k-l^ ds 



+ 



+ 



\ z \b ai { a Xl^Z k s+u )-b ai { a X k s -\^Z k s+ J i dr lai {u)ds 
\ z \b ai { a X k -\^Z k s+u ) - b aj ( a X k -\^Z k ^)\ d Va7 (u) ds 
J (g a (s, a X k ) g a ( 8 , "Xt 1 )) dWf 

■fff Vz\^{ a X k ^Z k +u )-p ai {<*X k -\^Z k +u ) 
1=1 Ja J-t 1 



7 =i- / -t 

I si 



£ f f Vzb ai { a X k -\^Z k s+u )-P ai {«X k -\^Z k -l) 

/=1 J0 J-T 



d Vaj (u)-dB^ 
dr) ai (u)-dB^ 



= f Af + Bf + C? + Df + Ef + F t a 



where we simply identify each of the six terms A t = (Af,a = 1, ■ ■ ■ , P), B u C t , D tl E t and F t with 
the corresponding expression in the previous formulation. By a simple convexity inequality (Holder) 
we have: 

\X k+1 - X k \ 2 < 6(\A t \ 2 + \B t \ 2 + \C t \ 2 + \D t \ 2 + \E t \ 2 + \F t \ 2 ) 

and treat each term separately. 

The term A t is easily controlled using Cauchy-Schwarz inequality, Fubini identity and standard 
inequalities and we obtain: 



to the Pd- 



E 


sup \A S \ 2 


< K 2 t [ E 


sup {Xt-Xt'l 2 


ds 




Lsu P s e[o,t] 


Jo 


L —T<U<S 





Similarly, the martingale term D t is bounded usir 
dimensional martingale (Jq (g a (s, a X k ) — g a (s, a ^S 



. P) and we obtain: 



E 


sup \D S \ 2 


< AK 2 


[ E 


sup \X k u -X k u ~'f 


ds 




L(Ks<t 




Jo 


1 T<U<S 





Let 



us now deal with the deterministic interaction terms B t and C t . We have: 



p f o 

/ (^z[b a7 ( a X k ^Z k +u )-b ai ( a X k -\"Z k +u )})d Vaj (u)d S 



i^i 2 = E 

a = l 7=1 

P ft P ,0 

^J2 tp / £/ ^z[\b a ,{ a X k s ^Z k +u )- 

a = l J 7 =l-/-r 



sup l^-x^- 1 ! 2 



< tP 2 L 2 f \X k - X k - % \ 2 ds < tP 2 L 2 f s 

JO Jo -t 

hence easily conclude that 

E[ sup \B S \ 2 } < tP 2 L 2 [ E[ sup \X k -Xl- l \ 2 ]ds 
se[o,t] Jo -t<u<s 
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and similarly, since Z and Z have the same law as X and X 



k-1 



*e[o,t] 



E[sup \C S \ Z ] < tP z L z I E[ sup {XZ-X^flds 



-T<U<S 



Eventually, the terms E t and F t use Burkholder-David-Gundy (BDG) inequality in place of Cauchy- 
Schwarz' together with similar arguments as used for B t and C t . Using the independence of the 
Brownian motions (_B" 7 ), BDG inequality yields, for the term E t (and similarly for the term Ft): 

E[ sup \O s \ 2 } < 4P 2 L 2 [ E[ sup {X^-X^-^ds 
se[o,t] jo -t<«<« 

for (9 t equal to _Ej or F t . Putting all these evaluations together, we get: 

el 



E 



sup \X 

se[o,t 

^k+l 



fc+1 



X 



k I 2 



<6{T + A)(K Z + 2P Z L z ) / E[ sup |X* - X* 

— T<U<S 



fe — 1 1 2- 



(4) 



Moreover, since X t +1 = X k for i G [— r, 0] by definition, we have, noting 



Af f = E 



sup \X k s +1 -X 

-T<S<t 



A' I 2 



the recursive inequality M t fc < if" f* M^' 1 ds with if" = 6(T + 4) (if 2 + 2P 2 L 2 ). We thus get by an 
immediate recursion: 



M t fe < (K") k 

(K»)*i' lf „ 



< 



o jo 



M° dsi . . . ds A 



(5) 



and M 4 ° is finite because of lemma [l] Routine methods starting from inequality ^ and allows proving 
existence and uniqueness of fixed point for <!> (see [211 PP- 376-377]). Indeed, the fact that 



sup \X k s +1 ~X k s \ 2 

se[-r,T] 



1/2 



< 00, 



implies the almost sure convergence of the series J^^Lo su Pse[-r,T] — X k \ an d hence the a.s. 

uniform convergence of the partial sums: 

n 

x? + Y,(Xt +1 *t) = x? 

k=0 

on [— r, T]. Denote by X t the thus defined limit. It is clearly almost surely continuous and JF t adapted. 
From lemma[l] we know that X 6 M 2 (C). It remains to show that X indeed satisfies the equation ([2|. 
This is easily done using the estimates we derived. Indeed, it is easy to show at this point that 



E[ sup \${X n ) t - $(X) t \ z ] < K" / E[ sup \X?-X a \']da->0 



-T<t<T 



-T<t<T 



as n — > oo. Hence <1>(X) is equal to the limit of the sequence <P(X n ), which is by definition, equal to 
X n+1 whose limit is X. This concludes the fact that X is a fixed point of ^ which completes the proof 
of the existence of solutions of the mean-field equation ^ . 
Uniqueness: 

Assume that X and Y are two solutions of the mean-field equations From lemma [l] we know that 
both solutions are in A4 2 (C). Using the bound Q, we directly obtain the inequality: 

rT 



E 



sup \X S -Y a \ 

T<S<T 



< K" 



E 



sup \X U -Y U \ 

-T<U<S 



ds 



which by Gronwall's theorem directly implies E 
a.s. follows. 



sup_ r<s<T \X S - Y s \ 2 = whence X = Y on [-r, T] 
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3 Convergence results 

We now turn to the main result, namely the convergence in law of the solutions of the network equations 
towards the equations More precisely, we are interested in showing that the propagation of chaos 
applies in this case. The propagation of chaos property states that provided that all neurons have 
a chaotic initial condition (i.e., i.i.d. in Ai 2 (C T )), then in the limit where the number of neurons is 
infinite, all neurons behave independently and the limit process satisfies the equation given by The 
proof follows standard proofs in the domain as generally done in particular by Tanaka or Sznitman 
28,26J. It is based on the very powerful coupling argument, that identifies the almost sure limit of 
the network solutions X t,N (the exponent N denotes the dependence of the law of neuron i upon the 
network size N) as the number of neuron tends to infinity, a method popularized by Sznitman in his 
extensive works (see e.g. [27]), the ideas of which traces back to the 70's (for instance, Dobrushin uses 
it in [H]). We start by considering the translation invariant case where the delay measures r)^ only 
depend on the neural population of neuron i, prior to dealing with the non translation invariant case. 



3.1 Quenched propagation of chaos for translation invariant neural fields 

We consider here rj^ = rj ai for all i in population a. We intend to prove the propagation of chaos 
property and the convergence towards the mean field equations for almost all realization of the the 
delays Ty under the assumption that for any fixed i, the random variables (Tj :) ')j g {i...jv} are independent. 
This is termed the quenched propagation of chaos property. In order to demonstrate this property, we 
define a coupling between the solutions of the network equations ([T]) and the mean-field equations ([2]) . 
Let i £ N such that p(i) = a. We define the process X 1 solution of the mean-field equation ([2]), driven 
by the Brownian motions (Wl) and (Bj 7 ) that govern the network process X 1 , and having the same 
initial condition as neuron i in the network, Q £ A4 2 (C T ). In details, X\ is the unique solution of the 
equation, for a = p(i): 

dX* t = f a (t, XI) dt + X£ =1 /° t M z [b ai {Xl Z] +s )]d Vai {s) dt 

+g a (t, XI) dWl + E^ =1 Jl° r Ez[/3« 7 (X|, ZJ +S )\ d Vaj (s) dB? t>0 , (6) 
XI =Ch{t) *G[-r,0] 

which constitutes a collection of independent stochastic processes (^t)i=i...jv that are identically dis- 
tributed with law X p ^ . Let us denote by mf the probability distribution of X" solution of the mean- 
field equation pi). As previously, the processes (Zj, • • • , Zf) constitute a collection of processes inde- 
pendent of (X|ji=i...jv and having the distribution m 1 ®- ■ -®m p . We aim at proving the almost sure con- 
vergence of a collection of processes (X^) := (Xl k ' , k = 1 . . . I) for some I £ N* and (h, - ■■ ,ii) £ N ; 
towards (Xl k ), implying the convergence of the law of (Xf ) towards N goes 

to infinity. We start by proving this result for / — 1 which will readily imply the case I > 1. In what 
follows, we shall denote by £ the expectation over the delays t^. 

Theorem 2 (Convergence in law) Under the assumptions (Hl)-(HJ^) and the chaotic initial con- 
dition assumption, the process [X l t ' ,—r < t < T), solution of the network equations ([T]), converge 
almost surely towards the process (Xj,— r < t < T) solution of the mean-field equations ([6]). This 
implies in particular convergence in law of the process (X t i,JV ,-t <t<T) towards (X?,-t < t < T) 
solution of the mean-field equations |2|. Moreover, if f and g are globally Lip schitz- continuous, we 
have for any i £ N and any T > 0: 



max £ { E 

i=l— JV 



sup \x^ n -x:< 2 



]} 



< °^ ( 7) 
min 7 (7V 7 ) v ; 



l -t<s<T 

where C(T) is a constant only depending on the parameters of the system and the time horizon T. 

Proof The proof differs from more classical proofs in that we consider that the network is composed 
of several distinct populations and includes random delays. Moreover, a special care has to be taken 
when dealing with the delayed interactions, and the fact that our drift and diffusion coefficients are 
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not globally Lipschitz-continuous, which again is not explicitly dealt with in the present proof but can 
be classically obtained using a localization argument. The principle of the proof consists in thoroughly 
analyzing the difference between the two processes as N tends to infinity. This difference is written as 
the sum of eight terms denoted A t (N) through H t (N): 



X l t ' N -Xl= I f a (a,Xy)-J.A,.X'J,l, / ,/,,(*. .VI- x )-//,. (.s..Y>/ir. 

p rt 



E/ w E (b« 7 (xi' N ,xi^-b aj (xixi'» Tij ))ds 
E fjr E (b a7 (xi,xi^)-b ai (xi,xi_ Tij ))d S 

7=1 J ° 7 p(j)=7 

E f ^ E b ai (Xl,XL TtJ ) - f MzK^Xi, Z] +a )} dr, ai {u) ds 



7=1 
P ,f 



P(j')=7 



7=1 J ° 7 p(j)=7 

E f w E (/3a 7 (x:,x^)-/3 Q7 (x:,xf_ Ti3 ))di?r 

7=1 1/0 7 p(j)=7 

E A^T E P^iXlXi )- f\ z [^(Xl,Z] +u )])d Va7 (u)dBl 



p(i)=7 

v4 t (A r ) + B t (N) + Ct(N) + D t (N) + E t (N) + F t (N) + G t (N) + H t (N) (8) 



Let us underline the fact that the probability distribution of these terms do not depend on the 
specific neuron i in population a considered. We are interested in the limit, as N goes to infinity, of 
the quantity £ {lEi[sup_ T<s< rp\Xl ,N — Xl\ 2 ]}. We decompose this expression into the sum of the eight 
terms involved in equation ([8]) using Holder's inequality and upperbound each term separately. The 
terms Af(N) and Bf(N) are treated exactly as in the proof of theorem [TJ and the control all the terms 
except E t and H t essentially uses on the same ingredients as in the proof of theorem [T] We denote K 
a common Lipschitz constant for / and g and deduce, similarly to the analysis performed in the proof 
of theorem [T] the following inequalities: 



E[sup < s < tAT[7 
E[sup < 

S<t/\TU 

E[sup 0< 
E[sup 0< 

S<tf\T V 

E[sup 0< 

S<tf\T V 

E[sup 0< 



\MN)\ 2 ] 

\B S {N)?\ 
\C S {N)\ 2 ] 
\Ds{N)\ 2 ] 

\G S (N)\ 2 ] 



R% T jtATv E[sup _ T<u<s lX i,N _ |2] du 
a V 2 r tAT u w r lr , I Yi.N vi |2l J„, 



^ J\ 1 J JC J [SUp_ T <„< s |yV„" 

<4X 2 / tA ^E[sup_ T<u< J^ 
< TL 2p^J^a nsup _ T<u<slx i, 

^ ^ L 2 P 2 f* ATU ■ - TC 1 (mm 



\X^-Xl\ 2 ]du 
X l u \ 2 }ds 



jo t 

< TL 2 P 2 J Q ~ ^[su P _ T <„< s |A - a;| 

'* Atu maXj = i...jvE[sup_ !; < u < g \X& N - Xl\ 2 } ds 

Maun , . I X i - N — X i \ 2 ] ds 



' • • * JO " .1- I.-"-* 

< 4L 2 P 2 /* E[sup_ T <„< s \X?" - Xl\'\ ds 

< AL 2 P 2 J* max i=1 ... w E[sup_ T < u < s \X^ N - Xl\ 2 ] ds 
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Let us for instance treat the case of F t (N) for the sake of clarity, all other terms are treated along the 
same lines. We have: 



E[ sup |^ s (iV)| 2 ] < PJ2 M 

7=1 



0<s<t 



sup 


I 


_0<s<t 





iv 7 



)dB. 



p(i)=7 

2" 



(BDG) < 4PE E 



7=1 
P ,t 



E ^ ( x s N > ^-i'-Tij ) _ p°<i (xi , x i'- Ti] ) 



PU')=7 



f 1 

(Cauchy-Schwartz) < 4PV] / — V 



E 



(assumption |(H2)P < AP 2 L 



7=1 J ° 
2 



p(j)=7 



/3 Q7 {X\' N , X 3 S '_ ) - /3 Q7 (X z s , X 3 S '" ) 



3,N 



ds 



< 4P 2 L 2 



E 



E 



X 



XI 



ds 



sup \X l u N - X 

-T<U<S 



i 2 



ds. 



In the decomposition {[sj) , the two terms E t (N) and -ffj (iV) are of a new form and remain to be eval- 
uated. Both term involves the difference between an empirical mean of a function of processes and an 
expectation term, and all have bounded second moment thanks to proposition [l] and assumption | (H3 ) [ 
We have: 

£{E[sup \E s (N)\ 2 ]} = £{m[suv I [ Y,W E (M^'^-ry)- / Vz[b an (Xi,Z2 +u )])d Vai {u)ds 

0<s<t <• 0<s<t Jo i iV 7 f~T J-r 

— — — — 7=1 0171=7 



<tpJ2 f £ {n^- E w*i.*f-T„)- fi* 

7=1 J ° 7 P0)=7 



7(^,^+ U ) d? ?a7( M )- 



(9) 



and using Burkholder-Davis-Gundy 

s{nsn V \H s (N)\ 2 ]} <apJ2 t £{m^- E 

S "* 7=1 70 7 P(J)=" 



/3 a7 (X],Xi_ Ti .)- / M z [(3 aj (Xl,Z] +u )dr, a7 (u)} 



ds 



Each of these expressions involve the processes Xf_ n . which, under the tensor product law of the Brow- 
nian motion and the delays, are independent and identically distributed population-wise. Moreover, 
we have: 



E 0{xi,xi- nj )- f®zie(xi,z] +u )}d Va ,(u)}\ 2 }} 

7 PU)=7 

1 ^ 

^Y, £ {^[( e ( X ^ X s-n 3 )-^Z,f a MXl 



N 2 
1 k,i=i 

where 6> € {b ail j3 ai }, noting that the expectation £{Ez[-]} of the composed random variable &{X\, X J S _ T .. ) 

under the law of X 3 and of the delays f = {t^} is precisely equal to J E^[6>(X], Z] +U )]dr) aj (u)]. 
In the above expression, f aj denotes a random variable with law r\ ai independent of the sequence of 
delays and of the Brownian motions. It is now easy to show that all the terms of the sum corresponding 
to indexes j and k such that the three conditions j 7^ i, k ^ i and j ^ k are satisfied are null. The 
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simpler way to see this property is to write the expectations as integrals with respect to the measure 
mf and observe that all terms annihilate, using the mutual independence of the processes (X^), of the 
random variables (Ty)j£{i...jv}> ancl the independence between the processes and the delays. Therefore, 
there are no more than 3 -/V 7 — 1 non-null terms in the sum (in the case a = 7 there are just iV 7 non-null 
terms), and moreover, all of these terms are uniformly bounded. The terms related to indexes j = k i 
satisfy the inequality: 



o{xi,xi)-M z ,,[e{xi,zi-)] 



}<2E{ 



E 



9(Xi,X>) + M Zt ,[0(Xi,Z]-)] 



<4K(1 + C 2 (s)) 



where C2(s) is the upperbound of the 1L 2 norm of X % given by lemmaJT] The terms related to the cases 
j = i (or symmetrically k = i) are bounded by the same constant, since we have for all k such that 
p(k) = a, by Cauchy-Schwartz inequality: 

sU [0(xi,xi_ Ti .)-M z , f ^[e(xi,z^j]Y ■ {0{xlx k s _ Ti .)-M z ^[0{xi,z^_ f j}) T I 



< Z{ E 



{' 



e{xixi_ Ti .)-v z , T [e{xiz«_ T ) 

<4^(l + C 2 (s)) 
We note C = 4RT(1 + C 2 (T)). We have: 



E 



0(Xi,X*_ Ti .) -M z , T [0(Xi,Z?_ T ) 



£<jE[ sup \E S (N)\- 

0<s<t 



> 7=1 



3AL - 2 , , 1 

7 <3T 2 p 2 C , 



N 2 

7 



min 7 (iV 7 ) 



and 



0<s<t 

Assembling all the estimates, we obtain 



£^E[sup \H S (N)\ 2 ] > < 12TP 2 C 



1 



min 7 (iV 7 ) 



max £{E[ sup |A^' JV - X l s \ 2 ]} < 8 < (K 2 + 2L 2 P 2 ) (T + 4) / max E sup \X' V - X'i 



i=l— jV 



0<s<t 



3=1— N l- T <u<s 



du 



3(T + 4)T— 



CP 



min 7 ( 7V 7 



Let us denote if' = 8(K 2 + 2L 2 P 2 ) (T + 4) and C = 24(T + 4)C P 2 . Since X\' N and X t ' are equal on 
[— r, 0], this inequality implies: 



max E[ sup \X\' N - X l s \ 2 ] < K' / max E sup \X 3 U - X^ 



i=l—N _ T < s < t 

Using Gronwall's inequality, we obtain 



j=l-N 



-T<U<S 



du 



a 



min 7 (iV 7 ) 



C 

max E[ sup \Xl - X\\ 2 ] < — — — . , 
i=i-..JV - T <s<t " mm 7 (iV 7 ) J 



z K ' C*-) ds < 



C 



K T 



K' min 7 (iV 7 ) 



(10) 



which tends to zeros as N goes to infinity by assumption. This property implies the convergence in 



law of (Xlf Tu ,-T <t<T) towards (X l t , -t <t<T). 

Corollary 1 (Propagation of chaos) Let I 6 IN* and fix I neurons (ii,--- ,ii) € N*. Under the 
assumptions of theoremlty the law of (XI 1 ' , • • • , X\ u , — r < t < T) converges towards m p ^ il ' ® ■ ■ ■ <S> 



m 
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Proof We have, using the notations of the proof of theorem [2] 



£ ^E 



sup (X t 

-T<t<T 1 



ii,N 



,X t 



ii,N 



i 

^E 

k= 

< 



E 



sup 




2' 


~T<t<T /\T V 







ic 



„K'T 



which tends to zero as N goes to infinity, hence the law of (X t 



ii,N 



X\ 1 ' , — r < t < T) converges 



towards that of {XI 1 , ■ ■ ■ ,X\\ —t < t < T A tjj). This property implies the convergence in law of 



N 



(8) m 



p(m) 



-t <t <T) towards (Xj 1 ,- ■ ■ , X\ l , -r < t < T) whose law is equal to rn p(il 



These two results readily yields the announced results on the quenched propagation of chaos and 
convergence for almost all realization of the delays in the translation invariant case. Let us now deal 
with the averaged convergence of the process in the non translation invariant case. 

3.2 Averaged propagation of chaos in the non-translation invariant case 



We now consider the case of non-translation invariant neural fields. In that case, the distribution of 
delays depends on the properties of the particular neuron considered (see e.g. the case of neurons in a 
box treated in section [4]). Let us consider neuron i in population a. The delays Ty between neuron i 
and neuron j of population 7 has a distribution given by rjij and these different laws, averaged across 
all possible choices of location of neuron i, has the law given by r\ ai . 

In the previous section, we were able to take into account possible correlations between delays: 
indeed, the only property used was that for any fixed neuron, the sequence of delays (Ty)j e {i...jv} 
were independent, whatever possible correlations between ry and tm for k ^ i. We make the same 
assumption here. We denote by £i the expectation over all possible choices of locations for neuron i, 
i.e. all possible distributions 77^. Following the same lines as done in the quenched heterogeneity case 
for translation invariant systems, we show the following: 

Theorem 3 Under the assumptions (Hl)-(H^) and the chaotic initial condition assumption, the pro- 
cess (X % t * , — r < t < T), solution of the network equations ([T]) averaged across all possible values of the 
delays converge towards the process (Xj, —t < t < T) solution of the mean-field equations ([6]). This im- 
plies in particular convergence in law of the process (£i[Xl^}, -t <t<T) towards (Xf, -r <t<T) 
solution of the mean-field equations ([2|. Moreover, if f and g are globally Lip schitz- continuous, we 
have for any i £ N and any T > 0: 



max £ 

i=l—N 



E 



sup |£j 

-T<s<T 



Xl' N ] - xi 



< 



C(T) 



(11) 



where C(T) is a constant only depending on the parameters of the system and the time horizon T. 

Proof The proof proceeds as the one of theorem [2] by thoroughly controlling the distance between 
£i[Xl' N ] and X\ under the norm ||Z|| 2 = E[sup_ T<t<T |Z S | 2 ]. This distance is decomposed into the 



sum of eight terms analogous to the quenched translation invariant case, that are controlled exactly 
in the same manner as done in the proof of theorem [2j except two terms E' t (N) and H' t (N) analogous 
to E t (N) and H t (N). Regarding these terms, we have: 



£{E[sup |^(iV)| 2 ]} 

0<s<t 



< 



£{e[ sup 

L 0<s<t 

P t* r 



p 

E 

7=1 



1 



E 

p(j')=7 



^z[b ai {X\ 



E[ 



1 



E 

P(j)=7 



£ % {b ai {Xl,Xi_ )}- / E z [6 a7 (X 



di] a7 (u) ds\ 2 ]\ 
2 

}ds] 



(«)] 



(12) 
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and using Burkholder-Davis-Gundy 



7=1 70 7 p(i)=7 



(«)] 



In that case again, the expectation of the process £i{f3 ai (Xl,Xl_ Ti .)} under X^ and the delays is 

J_ Mz\Pa-y(Xl, Z] +u ) dr) a j(u)], and hence developing this expression, we are left with at most 3iV 7 
terms that are bounded, and we can hence conclude on the convergence of the averaged law towards 
the mean-field equations, and obtain the announced speed of convergence. 



Let us now go deeper into the dynamics of these limit equations. 



4 Dynamics of the firing-rate model 



In the previous sections, we showed the convergence of network equations ([I]) towards well-posed 
delayed McKean-Vlasov equations These equations form a system of implicit equations in the 
space of stochastic processes. Understanding the qualitative dynamics of the networks in the limit 
N — > oo is hence highly challenging in this general context. In order to go deeper into the analysis of 
the dynamics of the mean-field equations, we consider the dynamics of the classical firing-rate model. 



4.1 Mean-Field equations for the firing-rate model 



We now instantiate the dynamics of our neurons as firing-rate models. In this model the intrinsic 
dynamics f a (x,t) = —x/9 a + I a (t) is affine, the diffusion function g a (x,t) is constant (equal to A Q ) 
and in which the interactions only depend on a nonlinear transform of the membrane potential of 
presynaptic neurons: b ai (x,y) = J a -yS(y) and /3 Q7 = o~ a ~S(y). It is obvious that the assumptions of 
the above sections are satisfied by the firing-rate model. Therefore, when considering the delays r/iy 
only depending on p(i), we have propagation of chaos and almost sure (quenched) convergence towards 
the mean-field equations: 



l+I a (t) + ^2 4 7 f M x [S(X^ s )}d V ^(s)) dt 

7=1 J - T 

P ,0 

+ X a dW t a + Y / ^- ( %[S(I 

7 = 1 J ~ T 



7 +a )]dT, ai {s)dBr (13) 



and if the delay distributions r/i 7 depend on the choice of i, the same result hold in an averaged sense. 
In these cases, we can reduce the dynamics as follows: 



Theorem 4 In the firing-rate model, solutions of the mean-field equations ( 13 ) are exponentially 
fast attracted towards Gaussian solutions. When the initial conditions are Gaussian processes, the 
solution is Gaussian, with mean fi a and variance v a satisfying a well-posed system of delayed differential 
equations: 



l 



2 

= -a-V, 



Ma + ELi S- T fiv-rit + s), u 7 (t + s))dr] aj (s) + I a (t) 

P ( \ 2 ( 14 ) 



where f((i, v) = M[S(X)] for X - J\f(ji, v) 
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Proof By the variation of constant formula we have: 

X? = X^ 6 * + J* e-C-*)/^ (j Q ( s ) + Jen f ^x[S(X] +u )]d V ^(u)) ds 
t p r t i-O 



! JO J-T 



The initial condition term IJe _( ' 9 ° vanishes exponentially fast. The other terms include deterministic 
terms and stochastic integrals with respect to Brownian motions of deterministic functions, hence are 
Gaussian. All but the initial condition term are hence Gaussian and the solutions are exponentially 
fast attracted towards Gaussian solutions. Taking the mean and covariance function of this process 



readily yields equations ( 14 1 



In the firing-rate case, we hence have an important reduction of complexity. We now exploit this 
simpler form of the mean-field equations to analyze the dynamics of the network using the theory of 
delayed differential equations (see e.g. |16j ) and compare it to numerical simulations. We particularly 
focus on the importance of the noise in the stochastic weights and on the distribution of delays. 



4.2 Distributed delays in a one-population network 

In order to analyze the impact of the distribution of delays on the dynamics of the network, we simplify 
the model considering a one-population network (P = 1, we hence drop the population indexes since 
there is no ambiguity here). We further assume that fj, = 0, or more precisely consider centered 
sigmoid functions (5(0) = 0) and no input (I = 0). In order to further simplify the system, we consider 
S(x) = ert(gx) where erf(x) = ^= Jq e~ x I 2 . In that case, integration by parts yields: 

.9M v 



f(p,v) = erf( 



v/i + 9 2 v' 



In that simplified case, it is easy to see that a stationary solution is given by (fi*,v*) — (0,A 2 /2). 
Due to the quadratic form of the variance equation, the linearized equation on the variance is trivial 
x a — —j-x a . Therefore, the level of noise in the connections do not come into play in that particular 
case, and the stability of the fixed point only depends on the characteristic equation related to the 
main, which is given by the dispersion relationship: 

£ = -\ + J9 f e^dr,{s). 

^/27r(l + <? 2 A 2 /2) J- T 

The solutions of this equations are the characteristic exponents of the system, and relate to the stability 
of the fixed point considered: if all characteristic exponents have negative real part, the equilibrium 
is asymptotically exponentially stable, and there exists a characteristic exponent with strictly positive 
real part, the equilibrium is unstable. Changes in the sign of the real part of the characteristic exponent 
constitute bifurcations of the system. Let us start by analyzing the presence of pitchfork bifurcations 
that correspond to real characteristic exponents crossing the imaginary axis. Such bifurcations hence 
occur when the following relationship is satisfied: 

o- -- Jg 

+ ^/2tt(1 + g 2 \ 2 /2) 

and hence this bifurcation occurs independently of the choice of the delay. Turing- Hopf bifurcations 
occur when the system has a pair of complex conjugate characteristic exponents with non-zero imag- 
inary part crossing the imaginary axis. Such bifurcations hence occur when there exists to > such 
that f = iui, i.e.: 

1 , Jg 



V^(1+ 5 2 A72) ' 

The existence of such bifurcations exist depending on the type of delay distribution. It is important to 
note at this point that there is no such bifurcation when delays are equal to zero. 
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4-- 2.1 Dirac delta delays 

Assuming that the neuronal populations are highly concentrated so that delays can be considered 
constant (i.e. r\ = 5- T ). The dispersion relationship can be easily solve and one finds that Hopf 
bifurcations exist only for 

m >!, 

^(i+5 2 f ) 

i.e. for strong enough connectivity and small enough noise. Under that assumption is then easy to show 
that Turing-Hopf bifurcations arise when the parameters satisfy the relationship, when J < 0: 

r = 7r arctan(w#) 



with ui = . / 9 , 2 1. These result into oscillations of the solutions at a pulsation equal by oj. 

V 2w(l+ fl ^) 1 4 3 

Details of the calculations are left to the reader and follow the lines of the proof used for localized delays 
of the next section. Note that the value of the delay corresponding to the Hopf bifurcation is a function 
of A. In particular, no Hopf bifurcation is to be found for A > A* = ^/2(J 2 /(27r) — 1/.9 2 )- The curve of 
Hopf bifurcations obtained in the plane (A, r) and simulations the moment equations are displayed in 
figure Fig.[l] These simulations illustrate the fact that, depending on the noise and the delay chosen, 




Fig. 1 Dirac delta distributed delays, g = 1, J = —2, A* ~ 2.449. For A = 0.5, the Hopf bifurcation occurs for 
a value of r ~ 1.33 and for A = 1, the Hopf bifurcation occurs for r ~ 1.727. 
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(a) A = 0.5, t = 1 (b) A = 0.5, r = 1.5 (c) A = 1, r = 1.5 

Fig. 2 Simulation of a 3 000 neurons network with Dirac delta distributed delays corresponding to the cases 
presented in Fig. [T] blue: 30 arbitrarily chosen trajectories in the network, and cyan dashed line: empirical 
average. The dynamics observed for the mean-field equations is clearly recovered in the network simulations. 



the law of the mean-field equations has a mean that is either periodic or constant, yielding periodic or 
stationary laws of the mean-field equations. Since all neurons have the same probability distribution, 
in the periodic case all neurons will oscillate in phase. These results are valid for numerical simulations 
of the stochastic network as displayed in figure Fig. [2] 

4-2.2 Localized uniformly distributed delays 

We now consider instead of Dirac delta distributions of the delays that the neurons in each population 
have a dispersion of order S around a typical value r. One of the simplest distribution modeling this 
kind of dispersion is the uniform distribution r\ — ~k[ T —s/2,T+s/2]dt/5. The Dirac delta delay can be 
seen as the limit of this distribution as 5 — s- 0. It is easy to show that we have in that case a dispersion 
relationship: 

1 Jg c _ gT e^/ 2 -e-^/ 2 

+ v^(l + <? 2 A 2 /2) e ' 
and possible Hopf bifurcations occurring when there exists uj > such that: 

iu) = -- + 2 Jg c -w sin M/ 2 ) 

^2tt{1 + g 2 \ 2 /2) 

Equating the real parts and the imaginary parts on both sides yields the system: 

= Jg sin^) cog(tjr) 

^/27r(l+3 2 A 2 /2) « V ' 

Jq sin(l7) . / \ 

y/2ir(l+g 2 \2/2) n K ' 

where we denoted Q = ujS/2. These equations are relatively hard to solve. But since we are interested 
in parameters (r, 5) corresponding to Hopf bifurcations, we can consider that these equations are self- 
consistent equations in the plane (r, 6) parametrized by the variable uj (or more precisely f2) and obtain 
the following equations on the Hopf bifurcation curves when J < 0: 



1 / 1 JV sin 2 (r2) 



£ 

" ~ ' 4f2 2 \ 2 ' 27r(l+.g 2 A 2 /2) Q 2 

7r — arctan(w6>) it — arctan(2J?^ / '6) 
uj ~ 2Q/5 

These curves are plotted for different choices of parameters in figure Fig. [3] For small absolute values 
of the connectivity coefficient J, the system presents a Hopf bifurcation curve which is an increasing 
function of 5: the larger the dispersion, the larger the necessary delay to trigger oscillations. For instance 
fixing J = — 2 and 6 = 1 and a fixed value of the delay r = 1.5, we show that increasing the dispersion 
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0.4 0.6 OS ) 12 . » „ „ „ „ s «„ 

(a) Bifurcation diagram (b) Bifurcation diagram 
J = —2, 8 = 1 J = -5, 6 = 2 




(c) r = 1.5, 5 = 0.6 (d) t = 1.5, S = 0.6, network 




(e) r = 1.5, 5 = 2 (f) t = 1.5, 5 = 2, network 

Fig. 3 Uniformly distributed delays in [r — 5/2, r + 5/2] for J = — 2 and # = 1 (all but (b)) or J = — 5 and 
9 — 2 (b), (? = 1 and A = 1. (a-b): Hopf bifurcations as a function of the spread of the distribution and the 
delay. Strange phenomena with multiple Hopf bifurcations can arise as J is decreased (b). (c-f): for r = 1.5, 
oscillations present in the Dirac-delta delayed network persist for small values of S. However, when S exceeds 
the value related to the Hopf bifurcation, the oscillations disappear. 

of the delays S destroys the oscillatory pattern: the network shows a transition from perfectly periodic 
in law synchronized oscillations to stationary asynchronous behaviors. For a different set parameters, 
we display a relatively complex bifurcation diagram presenting several parameter regions corresponding 
to the first Hopf bifurcation and hence to the presence of oscillations. 



4-2.3 Neurons in an interval 

Let us now further specify the type of distribution of delays that arises in neurosciences. As stated, 
delays are chiefly related to the axonal propagation at finite speed and to the specific time the signal 
transmission takes, t s . Assuming constant speed c, the delay r,-j between neuron i at location and 

a neuron j at location rj is = ^ r '~ rj ■ + t s . Delays are hence strongly related to the distribution of 
neurons in the cortex. Moreover, in that case, the system may not be translation invariant (if we do 
not assume for instance periodicity) as we will see below. We hence need to identify the distribution 
of distances between distributed points in specific topologies. Such distributions are known in closed 
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form for some simple topologies, such as the rectangle or the ball. But even in these favorable cases, 
the bifurcations appear harder to characterize. 

This is however possible rigorously in a 1-dimensional case. Assuming uniform distribution of the 
neurons in an interval [0,a] and that c = 1 (without loss of generality), it is easy to show that the 
distribution of the distance r between neuron i at location and other neurons has the density: 



, x dr 

Vi = (l[0,a-r ( ] + l[0,n 



a 



and the averaged law the density: 



'2 2r N 
dfj(r) = ( n)dr. 



a a 2 



In the case where we consider a periodic neural area (i.e. identifying and a), the distribution of 
distances is clearly the uniform distribution on [0, a/2] and the results of the previous section apply. 

When considering non-periodic neural areas, applying the results of the above section ensures that 
an averaged propagation of chaos takes place, and that the network equations converge towards the 
mean-field equations given by theoremUl A Gaussian process with mean and covariance A 2 /2e~^~ s ^ e 
is a stationary solution, i.e. the point (0, A 2 /2) is a fixed point of the moment equations ( [14] ). At this 
point the dispersion relationship now reads, denoting J = 



Jg 



y/2ir(l+g*\*/2)' 



Possible Hopf bifurcations occur when one can find u> £ R such that: 

1 .2J / .( \ e- iu 

iw = ---i— l+i — + 

a cua \ \Lua tea 

In order to solve this equation, we note [2 — cua and Z(f2) = ^ ^1 + i (jj + § -jj-j\ Remark that 

Z{Q) tends to 1/2 when a — > hence recovering the Dirac delta delay case. Equating modulus and 
argument in the Hopf dispersion relationship one obtains: 



a 2 = 



\Z{QW 



Ta =(~ + Arg(Z{f2)) - tan (^-j + 2fcvrj ^. 

Considering as above J7 as a parameter, we can find the locus of the Hopf bifurcations in the plane 
(a, t s ). We obtain a sequence of Hopf bifurcations indexed by k, and the first bifurcation is responsible 
for oscillations appearing in the system. Interestingly, we observe that the value of t s corresponding to 
oscillations takes a minimal value for a certain value of a (see figure Fig. [4). This can be understood 
heuristically as follows. In the interval [0, a], the mean distance between neurons is equal to | and 

2 

the standard deviation to Increasing a hence has the effect of linearly increasing the averaged 
effective delay, and quadratically the dispersion. For small values of a, the effective delay created by 
the dispersion of the neuron increases faster than the dispersion, and therefore a smaller constant delay 
t s is necessary to trigger oscillations. As a is further increased, the dispersion grows faster than the 
mean effective delay, and in that case the dispersion implies that larger delays are necessary to trigger 
oscillations, as seen in the uniformly distributed delays case. 

Fixing all parameters and changing the size of the interval in which the neurons lie can hence 
induces transition from stationary to periodic and back to stationary. Neurons are hence perfectly 
synchronized in a specific region of the neural field spatial extension, and otherwise asynchronous. 
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12 3 4 5 6 

(a) Bifurcation diagram 



(b) a = 0.1 (c) a = 1.5 (d) a = 3.5 

Fig. 4 Neurons uniformly distributed in [0, a] inducing propagation delay and with transmission delay t 3 : (a) 
Hopf bifurcation diagram in the plane (a,r s ) (same parameters as above). For r s = 1.1 (b-d), increasing the 
parameter a (the size of the neural field) induces transitions from stationary to periodic back to stationary. 



5 Conclusion 

In this article, we analyzed the limits and dynamics of large-scale neuronal network models in the 
presence of noise. We particularly focused on the presence of heterogeneous delays in the interactions 
between neurons, an ubiquitous phenomenon in neuronal network activity that is understood to have 
important effects on the brain function. We considered two different kinds of systems: (i) translation 
invariant systems where the delays have the same distributions whatever the neurons considered, 
and (ii) non translation invariant case where delays distribution depend on the particular neuron i 
considered. In both cases, we derived the mean-field equations associated with such randomly delayed 
networks for general models, and showed that the propagation of chaos occurred, in a quenched sense 
for translation invariant systems (i.e. for almost all realization of the delays), or after averaging across 
possible delay distributions in the non translation invariant case. 

In both cases, the limit equations are identical and given by a system of delayed McKean-Vlasov 
equations with distributed delays. The equations obtained are not classical stochastic equations: they 
involve a term depending on the law of the solution at previous times (due to the presence of distributed 
delays). We showed that these equations were well-posed, and admit a unique solution. Nonetheless, 
this solution remains relatively hard to characterize at that level of generality. In order to characterize 
the dynamics of the network and in particular the importance of heterogeneity in the interconnection 
delays, we instantiated a popular model in neuroscience, the firing-rate model. In the mean-field limit, 
the solutions converge towards Gaussian processes, whose mean and variance satisfy a closed set of 
differential equations with distributed delays. 

This property allowed, using the theory of bifurcations of delayed differential equations (see e.g. [Tt)]). 
to analyze in detail the role of delays in the macroscopic dynamics of neuronal networks. We showed 
that depending on the averaged delays, the network can either present a stationary or a synchronized 
periodic behavior. Such transitions as a function of the amplitude of delays were already observed in 
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heuristic models including constant delays. But beyond this fact, we showed that averaged delay is not 
enough to predict the behavior of the system, and a knowledge of the full delay distribution is neces- 
sary. We illustrated this fact exhibiting the bifurcation diagram of the Hopf bifurcations as a function 
of the variance of uniformly distributed delays, and probably more interestingly from the biological 
viewpoint the dependence of the network behavior as a function of the delays arising from the distance 
between the cells. We showed that the macroscopic behavior arising from such networks depend on the 
size of the neural area considered, that parametrizes the delay distributions in the network. This was 
done analytically in a simple one dimensional case, and observed that the size of the cortical column 
considered has non-trivial effects on the macroscopic dynamics of the neural area. 

This observation is very interesting from the neuroscience viewpoint. Indeed, it has been observed 
that across different species, the cortex either organizes into columns or shows a dispersed distribution. 
This is for instance the case of the primary visual areas: mice and rats for instance present no specific 
organization of the neurons responsible for orientation selectivity (a so called salt-and-pepper struc- 
ture), whereas the cat and primate present structured orientation selectivity columns (see e.g. [3T]). 
The reason of that difference is still poorly understood, and in that article, the authors suggest that 
the lack of orientation maps in rodent is related to the small size of VI in these species. They argue 
that making a synaptic connection is easier locally. This possibility of creating synaptic connections 
is not a definitive answer to that question since the authors report a few counter-examples such as 
the ferret or the tree shrew. Including the delay distribution in this picture and how the topology and 
size of the cortex influence the synchronization property at the level of cortical columns could allow 
going deeper into the analysis of the reasons of the functional reliability of orientation selectivity across 
different species with sensitively different cortex size. 

This work brings another refinement into the theoretical implication of delays in the dynamics of 
the network. The equations derived here provide a model of large-scale neuronal network with hetero- 
geneous delay, and a preliminary study of the dynamics of such networks was initiated. Enriching this 
model by considering several populations is a straightforward extension and the present manuscript, 
and analyzing such models would allow going deeper into the analysis of neuronal networks and on the 
effect of heterogeneous delays in relationship with specific cortical functions. Eventually, the cortex is 
far from being the only system presenting heterogeneous distributions of interconnection delays. Most 
networks, including communication networks, internet, social networks and physical networks such as 
spin glass present delays in the communications, that depend on the specific distance between nodes. 
This study can hence be applied to these other cases, and applications of that framework to the TCP 
protocol is currently under investigation. 
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